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Numerical Methods in Aero-Optics 


A Introduction 

This is the final report on the 20 month project “Numerical Methods in Aero-Optics”. The 
termination of funding for this project concludes 16 years of research under sponsorship of the 
Air Force Office of Scientific Research. A proposal for continuation of the work is pending with 
the U.S. Army Research Office. However, the principal investigator would also like to have the 
AFOSR consider continuation of our funding, relative to our record of interaction with Air Force 
Research Laboratories. 

Applied as an image of an object is formed, adaptive optics techniques compensate for degra¬ 
dations added along the path of the light from the object being imaged. Image restoration (post¬ 
processing) )tools are then used to scrub the captured optical image even cleaner. The first phase 
is a massive control problem. The second is a delicate inverse problem. Both adaptive optics and 
image restoration demand sophisticated mathematics and state-of-the-art computation. Figure 1 
gives an overall diagram of a typical adaptive-optics closed-loop control system. 



Figure 1: Typical Adaptive Optics System 

Image restoration involves the removal or minimization of degradation (blur, clutter, noise, 
etc.) in an image using a priori knowledge about the degradation phenomena. Blind restoration 
is the process of estimating both the true image and the blurring operator from the degraded 
image characteristics, using only partial information about degradation sources and the imaging 
system. Our main interest concerned optical image enhancement, where the degradation involves 
a convolution process. Image restoration techniques are also providing clearer views of objects. 
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The power of these tools is substantial. One of our simulations, for example, shows them 
improving the resolution of a telescope from being barely able to spot an object the size of a house 
trailer in earth’s orbit to detecting a hand waving from the trailer’s window! 

Adaptive optics compensation plays an essential role in current state-of-the-art atmospheric 
telescope imaging technology. The ideal earth-based astronomical telescope is built on bedrock, 
high on a remote mountain. The solid foundation partially stabilizes the telescope against wind 
and other potential causes of vibration, while the altitude and isolation minimize atmospheric 
degradation. The Hubble space telescope carries this logic to its natural extreme, but even the 
Hubble’s accuracy is limited by the effects of thermal stresses and other forces that shift the phase 
of the incoming light ever so slightly. 

Ideally, light from a distant object high above the earth’s atmosphere arrives at a telescope’s 
mirror as a single planar wavefront. The only limit on resolution should be diffraction by the 
telescope mirror aperture. In imaging through the atmosphere, tiny local variations in the index 
of refraction of the atmosphere induce small phase errors that make the incoming plane wave 
look more like a sheet of crumpled paper. The mirror then adds phase errors of its own; even 
a theoretically perfect min\« will be distorted by thermal stresses, not to mention the effects of 
small vibrations in the telescope structure. 

In this setting, active and adaptive optics attempt to compensate for these phase errors using 
as a reference the phase error in the image of a guide star, either a bright natural star very 
near the target image or a “star” created by directing a laser into the atmosphere. Guide stars 
are especially effective against the degradation by atmospheric turbulence of images collected 
by ground-based telescopes because they provide an estimate of the unknown blurring operator. 
Furthermore, thermal distortion and gravity can induce small deformations in lightweight mirrors. 
Active optics corrects these very low frequency errors by delicately nudging the primary mirror 
with hydraulic actuators. 

Adaptive optics corrects the higher frequency errors caused by atmospheric irregularities and 
telescope vibration. The distortion measured using the guide star drives a control system that 
adjusts a separate set of mirrors. (The primary mirror is too big to respond fast enough.) With 
adaptive optics, instruments like the 3.5-m telescope at the Starfire Optical Range of the U.S. Air 
Force Research Laboratory in New Mexico, can partially correct the image before it is recorded. 
Note that this real-time control requires extraordinarily high-speed computation - up to 10 billion 
floating point operations per second. 

Further, postprocessing restores the adaptive optics recorded image to a state even closer to 
perfection by filtering out any remaining noise and blur that can be distinguished from the image. 
The classic tool is regularized least squares; one of the newest techniques we investigated is based 
on the solution of a nonlinear partial differential equation. Like adaptive optics, both demand 
cutting-edge computation. 

In the following sections information is provided concerning our activities under this grant on 
these topics over the past 20 months. We believe that this work in computational mathematics 
has resulted in several contributions to the state-of-the-art in investigations into algorithms and 
software for high-performance, high-resolution imaging applications. Our work was particularly 
directed toward the development and implementation of innovative algorithms and techniques 
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for optical imaging. The following research topics (each involving the formulation and numerical 
solution of difficult ill-posed inverse problems) are described in the report: 

1. Adaptive Optics (AO): Atmospheric turbulence has been a limiting factor for imaging 
since telescopes were invented. Both atmospheric and telescope induced aberrations distort 
the spherical wavefront of arriving light. Adaptive optics is a scientific and engineering 
discipline whereby a distorted optical wavefront is compensated to correct for errors induced 
by the environment, e.g., a turbulent atmosphere, through which it passes. Our work in this 
area, joint with Brent Ellerbroek (DoD Starfire Optical Range), Paul Pauca and Xiaobai Sun 
(Duke), and Moody Chu (NC State), has made contributions to deformable mirror control 
algorithms, performance analysis, software, and DSP hardware design of on-line adaptive- 
optics systems. Additional applications include methods for focusing and propagating laser 
beams. 

2. Image Restoration (IR): We have tried to make novel contributions to the image postpro¬ 
cessing steps of restou tic it, enhancement and analysis of images taken through a turbulent 
medium such as the atmosphere. This work, joint with GRA Paul Pauca (Duke), Tony 
Chan (UCLA) and Curt Vogel (MSU), has applications in defense, including remote sens¬ 
ing, target identification and laser weapons programs, and to civilian technology, including 
astronomical, industrial, and medical imaging. We are concerned with fast nonlinear opti¬ 
mization methods for blind deconvolution. We considered novel algorithms for single and 
multiple imaging channels, based on phase and sensor diversity. 

Our purpose was to investigate the areas of adaptive closed-loop deformable mirror control 
systems, image postprocessing using blind deconvolution techniques and phase retrieval, and ac¬ 
curate corrections to the phase aberration problems encountered in radar systems. Parallelizations 
of the computational algorithms were developed and implemented on advanced architectures in 
conjunction with the DoD High Performance Computing Modernization Program. A description 
of some of our recent work, including color images, graphics and multi-media animations, is on 
our Web page at: 

http: / / www.wfu.edu / ~plemmons / 

under Research Projects. 
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B Objectives 

The objectives of this project were to conduct rigorous mathematical research on inverse problems 
arising in the areas of: (1) adaptive optics deformable mirror closed-loop control systems, (2) image 
postprocessing using blind deconvolution techniques and phase retrieval, as well as (3) accurate 
corrections to the phase aberration problems encountered in radar systems. High-resolution images 
are essential in many important applications in defense, science, engineering, law enforcement, and 
medicine. The need to extract meaningful information from degraded images is especially vital 
for such DoD applications as aero-optics imaging, surveillance photography, and modern synthetic 
aperture imaging systems. Sources of image degradation vary among application areas, but include 
atmospheric turbulence, turbidity in a fluid medium, motion blur, insufficient sampling, electronic 
noise, and numerous other effects. 

The goals of this project are summarized as follows: 

• Develop and evaluate innovative adaptive optics deformable mirror control algorithms, with 
concentration on performance analysis, software, and DSP hardware design of on-line adaptive- 
optics systems for imaging through turbulence. 

• Develop new iterative algorithms with effective regularization strategies for blind deconvo¬ 
lution optical image restoration and enhancement using adaptive filtering and multichannel 
phase diversity methods. 

• Exploit the general mathematical structure of computational problems found in several opti¬ 
cal imaging and other signal processing problems to develop a family of methods and routines 
that are adaptable to several scenarios in obtaining high-resolution images. 

• Develop a high performance software package, entitled “Parallel Image Processing Environ¬ 
ment (pipe)” , which utilizes the optical image processing algorithms developed in the course 
of this research. 

As we move toward the new millennium, the research results produced under this grant may 
very well have important impacts on science and engineering as part of a continuing development 
of the computational foundations of aero-optics technology. Some promising results and new ideas 
have been put forward and they indicate considerable potential for further progress in solving these 
important imaging problems in an efficient and stable way on modern computer architectures. The 
project resulted in a variety of new technologies in the form of robust and efficient algorithms as well 
as their implementation. Parallelizations of the computational algorithms were investigated and 
implemented on advanced architectures in conjunction with the DoD High Performance Computing 
Modernization Program. These were extensively tested on applications from optical imaging. 
Packaging the results of our research into reliable software will further facilitate the effective and 
timely transfer of vital new knowledge to DoD research laboratories and to industry. 
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C Major Accomplishments and New Findings 

In the following sections information is provided to help the AFOSR determine whether the re¬ 
sults of this DoD research project are consistent with U.S. Air Force missions. We believe that 
the work described here in computational mathematics has resulted in several contributions to 
the state-of-the-art in investigations into algorithms and software for high-performance, high- 
resolution imaging applications. Our work is particularly directed toward the development and 
implementation of innovative algorithms and techniques for optical imaging, including adaptive 
optics and image postprocessing. 

C.l Adaptive Optics Imaging Systems 

This research is concerned with several related projects in optical imaging. The projects to be 
described in this section concern adaptive optics control studies in collaboration with Dr. Brent 
Ellerbroek at the DoD Phillips Laboratory Starfire Optical Range. Much of the work is summa¬ 
rized in our papers [8, 12, 38, 39]. 

We envision that related applications of our work can be identified at, for example, the DoD 
Intelligent Optics Laboratory in Adelphi, MD. One important project concerns the efficient com¬ 
putation of covariance matrices that arise in the modeling of atmospheric parameters in adaptive 
optics (AO) systems. In a more general framework, our work concerns the development and im¬ 
plementation of fast and accurate algorithms for integral discrete transforms with applications in 
adaptive optics control systems, see [39]. We are currently investigating the potential of our tech¬ 
niques for use, for example, in adaptive optics deformable mirror control. A formal mathematical 
framework which unifies much of this work is given in our paper [8]. 

Many modern astronomical telescopes are now built with deformable mirrors that can be ad¬ 
justed dynamically. Real-time control of the separate actuators of such a mirror can accommodate 
distinctly different sources of error, such as wind shake and time-varying atmospheric distortion. 
Each source has its own characteristic temporal frequency; those of wind shake, for example, are 
typically much higher than those of atmospheric turbulence. The key to adjusting the mirror ac¬ 
tuators on the fly is to choose a basis set of mirror deformations, known as mirror control modes, 
which best control each disturbance at a bandwidth matched to its characteristic frequency. In 
contrast, correcting all disturbances at a common bandwidth will either allow high frequency er¬ 
rors to sneak through if this bandwidth is set too low, or add unnecessary noise to the correction 
of low frequency errors if the bandwidth is too high. 

Brent Ellerbroek of the Starfire Optical Range, principal investigator Plemmons, and other 
co-authors, have developed a multiple bandwidth modal control strategy that can minimize mean 
squared phase error in the image across multiple sources of error [12, 41]. These techniques advance 
previous approaches by enabling the simultaneous optimization of both the mirror control modes 
and the associated control bandwidths without choosing in advance the basis set of control modes. 

The optical performance of the system at a particular bandwidth is characterized by a matrix. 
The optimal control modes could be determined by finding the one unitary transformation that 
comes closest to simultaneously diagonalizing all of these optical performance matrices. However, 
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the approximate simultaneous diagonalization of more than two matrices is not an easy task to 
formulate algorithmically, and any such scheme would be computationally expensive. Instead, 
we use a novel trace maximization approach based on a hill climbing scheme relative to pairs of 
matrices in order to minimize mean squared phase error. 

In particular, we have studied in publications [12, 41] a non-smooth optimization problem 
arising in adaptive optics, which involves the optimal real-time control of deformable mirrors 
designed to compensate for atmospheric turbulence and other image degradation factors, such 
as wind-induced telescope vibration. The surface shape of this mirror must change rapidly to 
correct for time-varying optical distortions caused by these sources of image degradation. One 
formulation of this problem yields a functional f(U ) = £” =1 maxj{(U T MjU)u} to be maximized 
over orthogonal matrices U, where U and a fixed collection of n x n symmetric matrices Mj. We 
consider the situation which can arise in practical applications where the matrices Mj are “nearly” 
pairwise commutative. Besides giving useful bounds, results for this case lead to a simple corollary 
providing a theoretical closcu-Lrm solution for globally maximizing / if the Mj are simultaneously 
diagonalizable. However, even here conventional optimization methods for maximizing / are not 
practical in this real-time environment. The general optimization problem is quite difficult and is 
approached using a heuristic Jacobi-like algorithm. Numerical tests using the algorithm indicated 
that the performance of adaptive optics systems, such as those of interest to the Air Force, can 
be improved by the use of our Jacobi-like algorithm. 

These schemes can show substantial improvement over single bandwidth control methods. 
They also lend themselves to parallel implementation and real-time computing, see [12, 41]. 

C.1.1 Adaptive Deformable Mirror Control 

Specially designed deformable mirrors operating in a closed-loop adaptive optics system can par¬ 
tially compensate for the effects of atmospheric turbulence. The systems detect the distortions 
using either a natural guide star (point) image or a guide star artificially generated from the back 
scatter of a laser generated beacon. To illustrate the basic idea, we sketch the main components 
in an AO system are in Figure 2. These include the deformable mirror (DM), the wave front 
sensor (WFS), and the actuator command computer. Light in a narrow spectral band passing 
through the atmosphere is modeled by a plane wave. When traveling through the atmosphere 
that does not have a uniform index of refraction, light waves are aberrated and no longer planar. 
In a closed-loop AO system, this aberrated light is first reflected from the DM. Some of this light 
is focused to form an image, and some is diverted to the WFS that measures the wave front phase 
deformation. The actuator command computer takes measurements from the WFS and map them 
into real time control commands for the DM. How this translation is done depends on the criterion 
selected. 

In general, wave front sensing is a key aspect of many optical techniques and systems such 
as optical shop testing, interferometry and imaging through random media such as the earth’s 
atmosphere. We have used the massively parallel IBM SP2 at the DoD’s Maui High Performance 
Computing Center (MHPCC) for simulation studies in order to better understand the character¬ 
istics of our algorithms. In a paper [12] which appeared in the J. Optical Soc. Amer. A, we have 
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Figure 2: A simplified closed-loop AO system with main components. 


helped to develop a theory with applications to closed-loop adaptive control methods to adjust 
the shape of the telescope deformable mirrors in real-time. A second paper [41] was presented at 
the SPIE Annual Meeting in San Diego, July 1998. Numerical tests using this algorithm indicate 
that in the presence of windshake jitter, the performance of a closed-loop adaptive optics system 
can be improved by the selection of distinct and independently optimized control bandwidths 
for separate modes of the wave front distortion profile. Related interaction on minimal variance 
optimal wavefront reconstructor updating methods was discussed at February 1998 and March 
1999 visits to the DoD Starfire Optical Range in New Mexico. Work with potential applications 
to optimal reconstructor approximation methods was the objective of these visits. 

C.1.2 Modeling and Adaptive Optics System Design 

AO has emerged as the technique of choice to mitigate the blurring caused by atmospheric tur¬ 
bulence in large aperture imaging systems that allow extremely dim objects to be observed [23]. 
Modeling and evaluating the expected performance of AO systems is essential to the proper design 
of AO systems for large telescopes such as Gemini, Keck, Subaru, and VLT [11]. The objective 
of AO system design is to select hardware parameters and control approaches that will optimize 
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system performance for the expected operating conditions subject to the resources available, and 
maximize the resulting payoff for the intended scientific applications. The large range of possible 
observing scenarios and AO system parameters entails numerous cases to be considered, and fast 
computational approaches for performance evaluation and optimization are highly desirable. 

Modeling AO system performance is an essential part of AO system design. Among other 
modeling approaches, the linear systems model framework developed by Ellerbroek [13] provides 
first-order performance estimates which account for the full range of fundamental adaptive optics 
error sources and their interacting effects. However, such system performance evaluation and 
optimization models require intensive computations based on the covariance matrices for the 
statistical relationship between the WFS gradient measurements which drive the AO control loop 
and the turbulence-induced phase distortions to be corrected. Ellerbroek has given a general 
element-wise formulation for each ijth entry of the covariance matrices in the form 

<*ij = ci [jT C 2 n {h)dh^ j j Wi{x.)wj{y!) C 2 n {h) (1) 

X Io° ^ [ Jo {l^ S ) Jo (t^ S ) ~ X ] ds dh dx dx '> 


where Jo is the zero-order Bessel function of the first kind, and the quantities y and y' are scalar 
functions of x, x' € R 2 , and h € R. See Ellerbroek[14] for notation and further details. Denote 
the two-parameter Hankel transform [10] of a function f(x ) by 

poo 

h(a,b)= f(s)J 0 (as)J 0 {bs)ds, (2) 

Jo 


where a, b and s are variable nonnegative real numbers. For the inner integral with respect to s 
in (1), 


a = 


2tt y 


b = 


27iy 


(3) 


L 0 ’ " U 

The two-parameter Hankel transform is invoked for each and every covariance matrix entry in 
the linear system model of Ellerbroek[13, 14]. Due to the fact that y and y' are functions of x, 
x', and h, the efficient and accurate evaluation of the two-parameter Hankel transform becomes 
critical to the evaluation of the outer integrals in (1). For covariance computations, it is desired 
to efficiently evaluate (2) in a block-wise manner for N values of a and b, where N is proportional 
to a subset of deformable mirror actuators and wavefront sensor measurements, thus obtaining an 
N x N matrix 

H( a,b), a,b e R N . (4) 


Products involving the matrix H (a, b) are core operations for the integral computation in (1) after 
appropriate discretization. Figure 3 illustrates a particular example of the matrix H( a, b) that 
arises in adaptive optics simulations [13, 48]. 

The computation of the covariance matrices is intensive in two aspects. First, the computation 
for each given parameter set invokes the task of generating each and every entry of the covariance 
matrix, which is followed by calculations using the generated matrix for performance evaluation. 
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Figure 3: The two-parameter Hankel transform matrix H( a, b). 


Secondly, such computation is carried out many times over a large sampling space of AO perfor¬ 
mance parameters. To get a glimpse of the parameter space, we list some of the parameters in 
groups according to the physical meaning: (1) observing scenario parameters such as wavelength, 
aperture diameter, and zenith angle; (2) assumptions on the atmosphere such as the wind profile, 
the Cl(h ) (or turbulence) profile, and the turbulence outer scale; (3) architecture components 
of wave front sensing such as the WFS subaperture geometry, beacon geometry, pupil conjugate 
range and noise level; (4) architecture components of deformable mirrors such as the actuator 
geometry and conjugate range, and (5) level of hardware imperfection or limitations. 

In terms of the computation complexity, the size of the covariance matrices is proportional 
to the number of deformable mirror actuators and wavefront sensor measurements. Typically, 
the matrices are of order from 100 to 5000 for high-order AO systems designed to compensate for 
turbulence on very large telescopes. Each matrix entry invokes the evaluation of multiple integrals 
as seen in (1). Present computational approaches to evaluating each matrix element may sample 
the integrand at up to 10 4 points, where the integrand itself may be computed as a single or double 
infinite sum [14]. The subsequent computations with the covariance matrices require numerous 
inversions and multiplications of large matrices [13] to obtain performance estimates such as 
the residual mean-square phase error and the associated optical transfer function (OTF) of the 
telescope. To cover a parameter sampling space of a reasonable size, there may be hundreds or 
thousands of such covariance matrices to generate and compute with using today’s computational 
practice for performance evaluation. 

We have proposed novel approaches for efficient computations that exploit (i) the mathematical 
relationship among the entries of the matrix for each parameter set and (ii) the common computa¬ 
tion procedures and intermediate results among different sets of parameters. In preliminary work 
in [38], we proposed the use of fast Hankel transform methods [26, 27] combined with a matrix 
representation of H( a, b) in factored form, leading to fast computation of covariance matrices of 
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the form (1). Such matrix factorization formulation aids not only in revealing complexity and 
approximation accuracy (such as errors introduced by truncation and discretization), but also in 
identifying operations for which fast algorithms are known. A low rank matrix factorization may 
allow for efficient computation of matrix-vector products involving H( a, b) in the computation of 
covariances. Furthermore, its compact form representation also allows for storage savings. 

We derived in [38] the following factored form of H(a,b): 

H(a,b) = H T (a,b)+E, (5) 


where 

H t ( a,b) = MD u C t D n CD v M t . (6) 

The factors D u and D v are diagonal matrices, and C is the discrete cosine transform matrix. The 
middle factor D N is also diagonal and holds the tabulated values of f(s). The factor M is the 
one-dimensional equivalent of the matrix involved in particle simulations [20], 


M tj = l 


i 2 

0 , 


o2 ’ 


j < i 
otherwise 


(7) 


Note that the term E corresponds to a k-th order end-point quadrature correction and takes only 
0(N ) to apply and store. A significant advantage of our approach is that the factorization in (5) 
is given in terms of fast operators, all of which can be applied to a vector efficiently in at most 
0(N log N) operations. Hence, computing a matrix-vector product involving H(a, b) can also be 
done in O(NlogN) operations. 

A related fast approach for subsequent computations involving H (a, b) is based on a power 
series expansion [14]. In matrix form, the expansions can be done at different levels of granularity. 
We illustrate our partition template for the expansion in Figure 4. In other words, one can 
partition the matrix into smaller low rank subblocks which can be stored and applied efficiently. 
Using the ideas behind the fast multipole algorithm, the matrix can be partitioned into blocks 
whose size doubles as the blocks get further away from the diagonal blocks. Thus, the total cost 
for matrix-vector multiplication can be reduced to O(NlogN) using the hierarchical partitioning 
and rank-revealing scheme just described. 

We will continue to extend the ideas presented here with the aim of producing robust and 
modular software for efficiently forming and applying covariances matrices required in AO system 
design. Our preliminary studies by Pauca, Ellerbroek, Plemmons and Sun in [38, 39] have shown 
that novel fast Hankel transform techniques can reduce the computational complexity of currently 
used techniques by Ellerbroek [13, 14] from 0(N 4 ) to 0(N 2 log N). 
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Figure 4: Template of low-rank subblocks of H{ a, b) with / given in [13]. The numbers represent 
the numerical ranks of the subblocks. 


C.2 Image Restoration 

This research topic is concerned with the development and implementation of computational 
methods for the solution of certain large and often ill-posed inverse problems [15, 42] in image 
reconstruction and restoration, [2, 29], where the solution may not depend continuously on the 
data. In general, applications abound in science, engineering and medicine, for identifying and 
measuring useful features from image data, such as military targets, terrain features, handwritten 
alphanumeric characters, or medical artifacts such as tumors. This work also includes technologies 
for information management, data fusion, and visualization of extracted information. 

The techniques of adaptive optics are helping ground-based sensors to capture better, but 
hardly perfect, astronomical images. More mundane devices like surveillance cameras are cursed 
from the start with gritty, low-resolution images. In both settings, image restoration techniques 
can help obtain a clearer picture by separating the image from the degradations. 

Classic restoration techniques often model the received image as the sum of noise and a blurring 
operator acting on the true image. (The blurring operator is a convolution with the point-spread 
function that characterizes the aberrations.) These techniques seek to recover the correct image 
by solving an ill-conditioned least squares problem, usually regularized through the addition of 
some smoothness requirement to be satisfied by the restored image. 

This linear problem is a formidable computation because of the size of the blurring operator 
- its dimension is the number of pixels in the image - and its inherent ill conditioning. Reg¬ 
ularization can render the restoration problem solvable in practice, by neutralizing some of the 
ill-conditioning, but it also removes sharp edges and similar distinguishing features. Once distinc¬ 
tive characteristics, for example, can become unrecognizable. 

Leonid Rudin of Cognitech, Inc., and Stanley Osher of UCLA, formerly of Cognitech, have 
developed a widely used alternative known as the total variation (TV) technique because it mini¬ 
mizes the total variation of the image instead of its second derivative [44, 45]. Rather than insisting 
on a completely smooth image, TV requires only that it have bounded variation, permitting sharp 
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edges but eliminating spurious oscillations. They came to their approach in part through methods 
used for tracking shock fronts in gas dynamics calculations, a setting that also seeks to preserve 
sharp boundaries without introducing extraneous detail. 

The problem at the heart of a TV image restoration problem is equivalent to a nonlinear 
partial differential equation (the Euler equation of the constrained minimum variation problem). 
The restored image solves this steady-state problem. In their original work, Rudin and Osher found 
that steady-state solution by iterating in time from an appropriate initial condition. More recently, 
various researchers, including Tony Chan of UCLA, Curt Vogel of Montana State University, and 
principal investigator Plemmons, have proposed preconditioned conjugate gradient methods for 
iterating to the solution [1, 7, 30, 31, 33, 49, 50]. 

The down side of TV-based image restoration is its computational expense. For example, 
Rudin and Osher’s time-stepping scheme for solving the TV differential equation can be slowed 
by stability restrictions that force it to take relatively small time steps. 

In this respect, principal investigator Plemmons and his colleagues James Nagy of Southern 
Methodist University, Paul Pauca of Duke University, and Todd Torgersen of Wake Forest have 
proposed a compromise: Use TV methods to sharpen the estimate of the blurring operator ob¬ 
tained from an auxiliary source like a guide star, then apply quicker linear restoration algorithms 
to the image (see [30, 31, 33]). 

Since the blurring operator is localized, TV techniques can be applied to it fairly cheaply, 
thereby gaining from the ability of minimum total variation to preserve sharp transitions without 
the expense of a complete TV restoration. Using the TV estimate of the blurring operator, linear 
restoration is then applied adaptively to subregions of the full image allowing subregions of the 
image to converge at a more natural rate, a technique the authors call the space-varying regular¬ 
ization (SVR) technique. This novel SVR method reduces the possibility of excessive smoothing 
and magnification of noise during the linear restoration phase. 

Plemmons and his colleagues solve the linear restoration subproblems iteratively using SVR. 
The challenge is to continue the iterations long enough to amplify the components associated with 
the image but not so long that the noise is amplified as well. By simultaneously monitoring the 
size of the image components and the residual, or equation error, they can choose a stopping 
criterion appropriate for the portion of the image under study. This approach stops iterations 
early for a region with little spatial variation, such as an image of empty sky, but lets them run 
longer for a busier region that includes, say, a piece of the satellite or star under observation. 

Simulations using the 3.5-m telescope at the Starfire Optical Range show impressive improve¬ 
ments: Combined with the multiple-bandwidth adaptive optics control of Ellerbroek, Plemmons, 
and others, use of SVR strengthens resolution by a factor of about 50. A telescope that can discern 
nothing smaller than 30 meters at a range of 1000 km using its optics alone finds its resolution 
improved to 20 centimeters when a combination of adaptive optics and SVR postprocessing is 
used. 

Many image restoration problems fall into the broad category of blind deconvolution because 
it is necessary to estimate both the true image as well as the blur from the degraded image using 
only partial information about the blurring operator. Blind deconvolution approaches include 
those developed by Julian Christou of Phillips Laboratory and by Plemmons and his colleagues 
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Michael Ng of the Australian National University and Sanzheng Qiao of McMaster University (see 
[35, 36]). These particular approaches couple constrained optimization with nonlinear conjugate 
gradient methods. A major part of the the principal investigator’s recent image postprocessing 
work has been concerned with developing new, effective blind deconvolution methods. 

C.2.1 Current Status of our Restoration Work 

Image restoration and enhancement is an essential aspect of image analysis [2, 3, 12, 29, 43]. It 
involves the removal or minimization of degradation (blur, clutter, noise, etc.) in an image using 
a priori knowledge about the degradation phenomena. Blind restoration is the process of esti¬ 
mating both the true image and the blurring operator from the degraded image characteristics, 
using only partial information about degradation sources and the imaging system [28]. Our main 
interest in image postprocessing is in optical image enhancement (restoring the phase profiles 
of electromagnetic waves passing through a turbulent medium), where the degradation often in¬ 
volves a convolution process e.g., [43]. Our related interests include work on electromagnetic wave 
propagation, absorption, and scattering. 

The image formation process is often modeled as a Fredholm integral equation of the first kind: 

/ OO roc 

/ h(x,y;s,t)f(s,t)dsdt + ri(x,y ) (8) 

-oo J —OO 

where g(x,y ) is the observed (degraded) image, f(x,y) is the true (original) image (unknown), 
and rj(x, y) is assumed to be additive, Gaussian noise. Here, h(x, y; s, t ) represents the blurring 
point spread function (PSF). Our particular application concerns image restoration computations. 
The interest is on the development of fast algorithms and their extension to iterative blind decon¬ 
volution, where the blurring operator h as well as the image / is to be estimated. 

If we let T-L denote the blurring operator, then the image restoration problem associated with 
(8) can be expressed as a linear operator equation 

g = Uf + r). (9) 

Let u and v denote two-dimensional variables. If % is a convolution operator, as is often the 
case in optical imaging [12, 43], then the operator acts uniformly (i.e., in a spatially invariant 
manner) on /. Here, (9) can be written as 

Hf(u) = f h(u- v)f(v)dv , (10) 

where f2 is the domain of /. The problem is to both deconvolve and denoise the recorded image 
during the reconstruction process, and we refer to this as the demising and deblurring problem. 
In optical imaging, the kernel h in (10) is called the convolution point spread function (PSF). The 
Fourier transform of h is called the optical transfer function (OTF). After discretization of (9), the 
spatial operator H defined by h in (10) is a matrix that we denote by H. Here, in the spatially 
invariant case, H is block Toeplitz with Toeplitz blocks. Thus the fast Fourier transform (FFT) 
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can be used in computations involving H, e.g., [4, 5, 47], with efficient implementation possible 
on high performance architectures [19]. 

A classical approach employed for solving (9) is penalized least squares, also called Tikhonov 
regularization in the inverse problems literature [15]. It requires minimization of an expression 

mf-g\\ 2 + aJ(f), (11) 

where || • || denotes the norm on L 2 (D), a is a positive (regularization) parameter. The functional 
J(f ) serves the purpose of stabilizing the least squares problem and penalizing certain undesirable 
artifacts like spurious oscillations in the computed /. Various choices of J(f) can be made, 
including IIS'/!! 2 , where S is some smoothing differential operator, or the identity. This model, 
based on the Euclidean norm, leads to fast linear methods for computing /, and is often the method 
of choice by practitioners [2, 29]. However the use of other norms, such as the C\ norm, lead to 
nonlinear minimization methods which sometimes result in superior enhancement of blocky, noisy 
images, but with added computational cost. Rudin, Osher and Fatemi [45] have suggested an image 
enhancement method based on solving a nonlinear PDE constrained minimization problem where 
the function being minimized is the Total Variation (TV) of the image / = f(x,y). Numerous 
other papers on various TV approaches for image denoising or restoration have been written in 
the past five years, including, [1, 7, 30, 44, 49, 50]. We have chosen to use the numerical approach 
developed by Vogel and Oman [49, 50] for solving the nonlinear TV optimization problem. The 
outer iteration is a fixed point method and the inner iteration involves a preconditioned conjugate 
gradient method [22], applied to a large scale system 

(H*H + aL 0 (f k ))f k+l = H*g, k = 0,1,... 

For deconvolution, H is block Toeplitz with Toeplitz blocks. Here, a and /3 are chosen to stablize 
the computations. The matrix aLp(f k ) is a 2-D nonconstant Laplacian with five bands, associated 
with an elliptic operator. Its spectrum can vary widely over the outer iterations, i.e., with the 
index k, causing numerical difficulties that must be considered [49, 50]. 

Our recent work in [30, 31] concerns a new space-varying regularization (SVR) approach to 
image restoration using multiscale methods, and associated techniques for accelerating the con¬ 
vergence of iterative image postprocessing computations. Denoising methods, including TV min¬ 
imization, followed by segmentation-based preconditioning methods for minimum residual con¬ 
jugate gradient iterations [22], are investigated in our work. Regularization is accomplished by 
segmenting the image into (smooth) segments and varying the preconditioners across the segments. 
The method appears to work especially well on images that are piecewise smooth. The algorithm 
has computational complexity of only 0(£n 2 logn) per iteration, where n 2 is the number of pixels 
in the image and £ is the number of segments used. Also, functional and data parallelization of 
the algorithm is possible. 

Our approach is especially attractive for restoring images with low signal to noise ratios, and 
magnification of noise is effectively suppressed in the iterations, leading to a numerically efficient 
and robust regularized iterative restoration algorithm. These ideas are being extended to a novel 
Krylov subspace approach, to enhance one of our approaches to blind deconvolution. We began 


14 


another approach in [35, 36], using a nonlinear inverse filtering technique. We are also studying 
a phase diversity blind image deconvolution method in [6]. A recent SIAM News article [9] 
discusses, in part, some of our prior contributions in high-resolution imaging, reported in papers 
[25, 30, 31, 32, 33, 35, 36]. The additional feature of estimating the PSF as well as deblurring the 
image in blind deconvolution adds another level of difficulty to the restoration process, leading to 
some exciting new challenges. 

C.2.2 Novel Blind Deconvolution Methods 

A fundamental issue in image restoration is blur removal in the presence of observation noise, using 
only partial information about degradation sources and the imaging system, i.e., the blurring 
operator is essential unknown. In the important case where the blurring operation is spatially 
invariant, then the basic restoration computation faces the usual difficulties associated with ill- 
conditioning in the presence of noise [3]. The image observed from a shift invariant linear blurring 
process, such as an optical system, is described by how the system blurs a point source of light 
into a larger image. The image of a point source is the point spread function , denoted by h. The 
observed image g is then the result of convolving the PSF h with the “true” image, say /. This 
blurring process is represented by the convolution equation g = h*f. The standard deconvolution 
problem is to recover the image / given the observed image g and the blurring operator h. The 
PSF of an imaging system can sometimes be described by a mathematical formula. More often, the 
PSF must be estimated empirically. Empirical estimates of the PSF can sometimes be obtained 
by imaging a relatively bright, isolated point source. In astro-imaging, the point source might be 
a natural guide star, or a guide star artificially generated using range-gated laser backscatter, e.g, 
[12, 43]. Notice here that the PSF as well as the image may be degraded by noise. 

Recursive Inverse Filtering In many applications data corresponding to h is not completely 
known. Blind deconvolution is the process of estimating both the true image / and the blur h 
from the degraded image g. One direction of our recent work on blind deconvolution, reported in 
[35, 36], is to incorporate regularization into and refine a nonlinear recursive inverse filter blind 
deconvolution method first proposed by Kundur and Hatzinakos [28]. They call their scheme 
the nonnegativity and support constrained, recursive inverse filtering method , or NAS-RIF, for 
short. (Also, some important direct filtering approaches to iterative blind image deconvolution 
have recently been proposed in [7, 16, 24, 53]). 

An algorithm for regularized iterative blind deconvolution using truncated eigenvalue and total 
variation regularization in conjunction with inverse filtering is developed in our paper [35]. The 
scheme is based upon a constrained optimization method, using a minimization procedure with 
nonlinear conjugate gradients (see [37], Chapter 5). We apply regularization to the inverse filter 
by using an inexpensive eigenvalue truncation scheme, and allow the user the option of applying 
total variation regularization to the estimated image. We call our approach the nonnegativity and 
support regularized recursive inverse filter (NSR-RIF) algorithm. The constrained optimization 
problem associated with our algorithm is convex-, thus, a variety of nonlinear optimization schemes 
can be used to minimize the objective function. See Table 1 for an outline of the proposed algorithm 
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and Figure 5 for an schematic diagram. We remark that the regularized image estimate can replace 
y as an input to compute the error vector z in the recursive inverse filter algorithm. The switch 
in Figure 5 indicates this option. 


Space Object Identification by Nonlinear RIF Blind Deconvolution 


f = g *s y*=f = “true” image g = “observed” image s = “inverse” PSF 



True Immge Restored Image 

[ Unknown f) 


Figure 5: Diagram of the NSR-RIF algorithm. 

Numerical tests are reported in [35, 36] on some simulated and optical imaging problems, and 
a comparison is made with the NAS-RIF algorithm [28]. Figure 6 also gives a sample restora¬ 
tion of a ground-based image of a satellite using our NSR-RIF blind deconvolution method from 
[35]. In this example a 256 x 256 image is considered. Specifically, the true object is an ocean 
reconnaissance satellite. A computer simulation algorithm was used to produce an image of the 
satellite, shown in Figure 6(left), as it would be observed from a ground-based telescope using 
adaptive-optics compensation. The satellite was modeled as being 12 meters in length and in an 
orbit 500 kilometers above the surface of the earth. The charge-coupled device (CCD) for form¬ 
ing the image was a 65,536 pixel square array. CCD root-mean-square read-out noise variance 
was fixed at 15 microns per pixel to reflect a realistic state-of-the-art detector. The computed 
restoration is shown in Figure 6(right). 

Recently there has been growing interest and progress in using total least squares (TLS) meth- 
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Definitions: 

g: the blurred and noisy signal of size k. 
n: the support size. 
p: the filter length is 2p + 1. 

s(t): the FIR filter parameter vector of dimension 2p + 1 at the t-th iteration. 

S(t): the corresponding matrix of the FIR filter s(t). 

y(t ): the estimate of the original signal at the t-th iteration. 

J re g(t): the objective function at s(t). 

XjJ reg {t) : the gradient vector of J reg (t). 
e: the tolerance for the termination. 

Initial Conditions: 

Set s(0) to all zeros with a unit spike in the middle, 
or, to an estimate corresponding to the inverse of the PSF. 

Iterations: 

For t = 0,1,... 

Compute y(t) = S(t)g. 

Project y(t) onto y(t) for nonnegativity. 

Compute the error z(t ) = y(t) — y(t). 

Compute s(t) = Fs(t), where F is the DFT operator. 

Project s(t) onto s(t) for finite support and TV regularization. 

Compute the error v(t) = s(t) — s(t). 

If J re g(t) < c, then stop; otherwise compute \/J reg (t)- 
If t = 0, then set the conjugate gradient direction vector d(t) = — V Jreg(t)', 
otherwise compute e(t) = [\/J reg (t) — VJreg(t — 1)]* V Jreg(t)/\\ V ^rej(^)||2) 
and set the direction vector d(t) = — V Jreg(t ) + e(t)d(t — 1). 

Perform a line minimization to determine S t such that 
Jreg(s{t) + Std(t )) < J reg (s(t ) + 8d(t )), \/5 G TZ. 

Compute s(t + 1) = s(t) + 8 t d(t). 

End For 


Table 1: Proposed NSR-RIF Inverse Filter Algorithm for Blind Deconvolution. 
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Figure 6: Degraded Image (left) and Restored Image (right). 


ods for solving blind deconvolution problems arising in image restoration. Here, the true image 
is to be estimated using only partial information about the blurring operator, or point spread 
function, which is subject to error and noise. In this respect, we have contributed a new iterative, 
regularized, and constrained TLS image restoration algorithm. Neumann boundary conditions are 
used to reduce the boundary artifacts that normally occur in restoration processes. Numerical tests 
are reported on some simulated optical imaging problems in order to illustrate the effectiveness of 
the approach, as well as the fast convergence of our iterative scheme, as indicated in our recent 
paper [34]. Preliminary numerical results indicate the effectiveness of the method. Future work 
on this project will exploit the constrained total least squares technique in phase diversity-based 
deconvolution arising from astronomy, extending work by Chan, Plemmons, and Vogel [6]. 

Our next steps are also to extend the phase diversity, total least squares and nonnegativity and 
support regularized recursive inverse filter blind deconvolution methods in [6, 34, 35, 36] to allow 
for multiple frames of data, where a large number of short exposure images are obtained over a 
short time period, as in video imaging [2]. This work overlaps with our research topic discussed 
next - multiframe (multichannel) blind deconvolution. 

Multichannel Blind Deconvolution and Phase Diversity Another approach to blind de- 
convolution that we have pursued capitalizes on temporal or sensor diversity to obtain stable 
reconstructions with enhanced resolution. It is based on a multichannel formulation, where noisy, 
differently blurred versions of the same image are available. In particular, sensor diversity arises in 
remote sensing, where the same scene is observed at successive time instants through the turbulent 
atmosphere, with a different transfer function at each instant (this can include the effect of partial 
compensation by adaptive optics [12, 43]). In some cases, the transfer functions differ primarily 
by their phase, giving rise to sb called phase diversity. Phase diversity involves the simultaneous 
collection of two (or more) short-exposure images. One of these images is the conventional image 
that has been blurred by unknown aberrations. An image is collected in a separate channel, by 
blurring the first image by a known amount, e.g., using a beam splitter and an out-of-focus lens, 
which is a quadratic blur. It is somewhat remarkable that estimates for the object and the un- 




known aberrations can be mad; from these two images. This is a research topic of considerable 
current interest in the imaging community, see [6, 9, 40, 46]. 

Using the two images collected by the phase diversity method as data, one can set up a math¬ 
ematical optimization problem (typically using a maximum likelihood formulation) for recovering 
the original image as well as the phase aberration. In the seminal work reported by the authors 
of [40, 46] an expectation-maximization algorithm is constructed to recover both the aberrations 
(blur) and the fine-resolution image common to both images. 

The phase diversity procedure is briefly described as follows [40, 43, 46]. The incoherent 
isoplantic image formation process is well modeled by the following convolution (see [43]): 

9jk(u) = Ujkiu - v)f(v) + T] jk , ( 12 ) 

where / is the object array, 'Hjk is an incoherent point spread function, corresponding to the jth 
atmospheric realization and the A;th diversity channel, gj k is the corresponding observed image, and 
u and v are 2-dimensional coordinates. The data set should be chosen to contain sufficient image 
frames (data) from a total of J atmospheric realizations. But, typically, only K = 2 diversity 
channels are needed. Phase diversity is introduced in the system’s coherent transfer function (see 
[40]), 

Pjt(z) = P(z)e |i ^M+».M 1>, (13) 

where P(z ) is a binary function that serves as a model of the telescope pupil function, 4>j{z) is 
the unknown phase-aberration (blurring) function associated with the jth. atmospheric realization, 
0 k (z) is a known phase diversity function associated with the fcth diversity channel, and z is the 
discrete spatial-frequency variable. The incoherent point spread function 'Hjk in (12) is just the 
squared modulus of the Fourier transform of the coherent transfer function given by (13). The 
objective is then to use (12) and (13) to set up an optimization problem for recovering the original 
image / as well as the phase aberration function 4 i. 

The expectation-maximization algorithm is often used to recover both the aberrations 4>j and 
the fine-resolution image / common to the images [46]. This is a computationally intensive al¬ 
gorithm. Our beginning work in [6] aims at improving the numerical efficiency of the restoration 
procedure. We have begun a project to develop fast computational algorithms for this phase 
diversity-based blind image deconvolution approach. In particular we addressed the following 
problems: 

• Regularization of the phase function in (13). This improves the stability of the mathematical 
optimization problem and eliminates some of the spurious minima. 

• Quasi-Newton methods, including Gauss-Newton, BFGS and nonlinear EN-like methods 
[51, 52], that take advantage of the special structure of the optimization problem. These 
combine fast convergence with relatively low cost per iteration. 

• Effective parallel implementation of multilevel preconditioners to quickly solve the large 
structured linear systems which arise at each quasi-Newton iteration. 
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• A theoretical and parallel performance analysis of the interaction of particular iterative 
methods for nonlinear and linear system solving, with the preconditioner and architecture 
e.g., combinations like nonlinear EN-like methods, linear EN-like methods, and various pre¬ 
conditioners [51]. 

Multichannel blind deconvolution techniques present both computational and theoretical chal¬ 
lenges. They can require solution of large least-squares problems for structured matrices that are 
block Toeplitz with Toeplitz blocks. Because of the huge sizes of some of the matrices arising in 
the problems (e.g., kn 2 x n 2 for k n x n frames of images, where 1024 x 1024 images are not uncom¬ 
mon), careful implementations of iterative algorithms that avoid forming these matrices explicitly, 
and which use preconditioning for speedup, are required. We developed methods incorporating 
statistical considerations to improve robustness to noise. Likewise, in the future, we expect to 
develop methods for model order reduction, to both improve algorithm stability, and to offer a 
tradeoff between performance and computation (see [17,18, 21]). We have also explored further the 
incorporation of nonlinear (and sparsity-based) regularization techniques [30, 31, 32, 33]. These 
steps provide yet additional ways to further enhance our methods for obtaining high-resolution 
aero-optics images. 
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Collaborative Research and Transitions at U.S. Air Force 
Laboratories 


Recent activities for-this work included visits to the Air Force Research Laboratory (AFRL) and 
Starfire Optica Range, Kirtland AFB, NM and to Wright Laboratory, Wright-Patterson AFB, OH. 
We are also in contact with DoD researchers at the Maui Space Surveillance Center (MSSC) in 
Hawaii. The primary contact at the AFRL is Dr. Brent Ellerbroek at the Starfire Optical Range, 
and we are collaborating on research involving closed-loop adaptive-optics systems. In two papers 
the authors have helped to develop a theory with possible applications for closed-loop adaptive 
control methods to adjust the shape of these mirrors in real-time [12, 41] . An paper on adaptive 
optics systems performance evaluation has been completed [39]. An abstract, “Leading Edge 
Methods in Optical Imaging”, was prepared for the DOD publication, Success Stories in High 
Performance Computing in 1998, concerning our use of the Maui High Performance Computing 
Center’s IBM SP2. Such collaboration is continuing, and a project involving optimal minimal 
variance estimators in adaptive-optics was just been completed in March 1999 [8]. 

We have also interacted with Dr. Julian Christou, Dr. David Tyler, and Dr. Donald Washburn 
at the AFRL on blind deconvolution and the airborne laser weapons program, and with researchers 
at Maui Space Surveillance Center (MSSC) concerning use of our methods for near real-time 
implementation in conjunction with their speckle imaging system. The MSSC has a direct line 
connection to the Maui SP2. Air Force Capt. Bruce Stribling is the primary contact at MSSC. 
Dr. Ellerbroek and Dr. Tyler at the AFRL have also supplied us with DoD satellite image data, 
sample phase screens, etc., for tests with our image post-processing work. Overall, we will be 
making available a robust, scalable software system that is portable across a variety of massively 
parallel computers and clusters of workstations. 

Technology transfer of our work in relation to the projects at DoD Laboratories has included 
research on ground-based imaging of satellites, and related aero-optics activities. The PI has par¬ 
ticipated in three DoD workshops at the AFRL in New Mexico: (1) the Aero-Optics and Image 
Reconstruction Workshop concerning airborne laser weapons development, (2) the Smart Sensors 
Workshop concerning remote sensing for wide angle satellite surveillance, where the satellite sys¬ 
tems will have on-board image processing capabilities when deployed, and (3) a second Workshop 
concerning recent trends in missile tracking for the airborne laser weapons program. Some of 
our work in these research projects concerns real-time adaptive filtering methods. Applications 
include closed-loop active noise (vibration) cancellation, with the potential for stabilizing firing 
pads for laser weapons. 

In addition to the DoD applications described above, potential technology transfer of our 
research on these imaging projects to civilian technology include astronomical and medical imaging, 
and remote sensing for commercial purposes. In astronomical imaging there is technology transfer 
to the astronomical community at large. Medical imaging applications of the adaptive-optics 
work include fluorescence microscopy in three dimensions, and the use of adaptive-optics methods 
as an aid in deblurring images of the retina through the eyeball. Our image post-processing 
methods can also be applied to enhancing satellite collected images of the earth for agricultural, 
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law enforcement, and geophysical purposes. 

During a visit by Plemmons to the DoD AFRL and Starfire Optical Range, February 25-27, 
1998, three presentations were given: (1) Fast Algorithms for Phase-Diversity Based Imaging, (2) 
Adaptive Condition Estimation for Image Reconstructor Updating, and (3) Image Postprocessing 
by Constrained Inverse Filtering. 

In addition, Plemmons and GRA Pauca, and others, visited the AFRL in New Mexico again 
on March 7-10, 1999. The purpose was to discuss the following topics in a workshop organized by 
Dr. Ellerbroek: 

• Control algorithms for real-time adaptive optics (AO), performance optimization for a priori 
conditions, and adaptive methods for unknown conditions. 

• Efficient numerical methods for AO system evaluation and optimization. 

• Early tests on full amplitude and phase conjugation in aero-optics imaging from aircraft. 

• Efficient algorithms for multiframe blind deconvolution and phase recovery for image restora¬ 
tion. 

This latest visit was taken to this Air Force Research Laboratory without AFOSR support , which 
ended 12/31/98. Our purpose in participating in the workshop was to facilitate the completion 
of some joint papers involving Dr. Brent Ellerbroek at the Starfire Optical Range, and to meet 
with other researchers on optical imaging, as part of an informal meeting. We have a proposal 
pending to the ARO, and we envision that related applications of our work can be identified at, 
for example, the Army Research Laboratory at Adelphi, MD. 

I Inventions or Patent Disclosures 

None. The research sponsored by the AFOSR under this grant concerned the development 
of rigorous mathematical models, computational algorithms and high performance computer soft¬ 
ware implementations that will be made readily available to the DoD and the private sector. 
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